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Abstract 

In this paper we apply an extended Kalman filter to track both the location and the 
orientation of a mobile target from multistatic response measurements. We also analyze 
the effect of the limited- view aspect on the stability and the efficiency of our tracking 
approach. Our algorithm is based on the use of the generalized polarization tensors, 
which can be reconstructed from the multistatic response measurements by solving a 
linear system. The system has the remarkable property that low order generalized 
polarization tensors are not affected by the error caused by the instability of higher 
orders in the presence of measurement noise. 
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1 Introduction 

With each domain and material parameter, an infinite number of tensors, called the Gen- 
eralized Polarization Tensors (GPTs), is associated. The concept of GPTs was introduced 
in [6, 4]. The GPTs contain significant information on the shape of the domain [3, 7, 9]. 
It occurs in several interesting contexts, in particular, in low- frequency scattering [15, 4], 
asymptotic models of dilute composites (see [23] and [10]), in invisibility cloaking in the 
quasi-static regime [8] and in potential theory related to certain questions arising in hydro- 
dynamics [24]. 

Another important use of this concept is for imaging diametrically small conductiv- 
ity inclusions from boundary or multistatic response measurements. Multistatic response 
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measurements are obtained using arrays of point source transmitters and receivers. This 
measurement configuration gives the so-called multistatic response matrix (MSR), which 
measures the change in potential field due to a conductivity inclusion. In fact, the GPTs 
are the basic building blocks for the asymptotic expansions of the perturbations of the MSR 
matrix due to the presence of small conductivity inclusions inside a conductor [17, 12, 6]. 
They can be reconstructed from the multi-static response (MSR) matrix by solving a linear 
system. The system has the remarkable property that low order generalized polarization 
tensors are not affected by the error caused by the instability of higher orders in the pres- 
ence of measurement noise. Based on the asymptotic expansion, efficient and direct (non- 
iterative) algorithms to determine the location and some geometric features of the inclusions 
were proposed. We refer to [4, 5] and the references therein for recent developments of this 
theory. An efficient numerical code for computing the GPTs is described in [11]. 

In [2], we have analyzed the stability and the resolving order of GPT in a circular 
full angle of view setting with coincident sources and receivers, and developed efficient 
algorithms for target identification from a dictionary by matching the contracted GPTs 
(CGPTs). The CGPTs are particular linear combinations of the GPTs (called harmonic 
combinations) and were first introduced in [8] . As a consequence, explicit relations between 
the CGPT of scaled, rotated and translated objects have been established in [2], which 
suggest strongly that the GPTs can also be used for tracking the location and the orientation 
of a mobile object. One should have in mind that, in real applications, one would like to 
localize the target and reconstruct its orientation directly from the MSR data without 
reconstructing the GPTs. 

In this paper we apply an extended Kalman filter to track both the location and the 
orientation of a mobile target directly from MSR measurements. 

The Extended Kalman Filter (EKF) is a generalization of the Kalman Filter (KF) 
to nonlinear dynamical systems. It is robust with respect to noise and computationally 
inexpensive, therefore is well suited for real-time applications such as tracking [26]. 

Target tracking is an important task in sonar and radar imaging, security technologies, 
autonomous vehicle, robotics, and bio-robotics, see, for instance, [13, 14, 16, 18, 19, 25]. An 
example in bio-robotics is the weakly electric fish which has the faculty to probe an exterior 
target with its electric dipole and multiple sensors distributed on the skin [1]. The fish 
usually swims around the target to acquire information. The use of Kalman-type filtering 
for target tracking is quite standard, see, for instance, [13, 14, 16, 18, 19, 25]. 

However, to the best of our knowledge, this is the first time where tracking of the 
orientation of a target is provided. Moreover, we analyze the ill-posed character of both the 
location and orientation tracking in the case of limited-view data. In practice, it is quite 
realistic to have the sources/receivers cover only a limited angle of view. In this case, the 
reconstruction of the GPTs becomes more ill-posed than in the full-view case. 

It is the aim of this paper to provide a fast algorithm for tracking both the location and 
the orientation of a mobile target, and precisely analyze the stability of the inverse problem 
in the limited-view setting. 

The paper is organized as follows. In section 2 we recall the conductivity problem and 
the linear system relating the CGPTs with the MSR data, and provide a stability result in 
the full angle of view setting. In section 3 we present a GPT-based location and orientation 
tracking algorithm using an extended Kalman filter and show the numerical results in the 
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full-view setting. In section 4 we analyze the stability of the CGPT-reconstruction in the 
limited-view setting and also test the performance of the tracking algorithm. The paper 
ends with a few concluding remarks. An appendix is for a brief review of the extended 
Kalman filter. 

2 Conductivity problem and reconstruction of CGPTs 

We consider the two-dimensional conductivity problem. Let B be a bounded C 2 -domain of 
characteristic size of order 1 and centered at the origin. Then D = z + 5B is an inclusion 
of characteristic size of order 5 and centered at z. We denote by < k ^ 1 < +oo its 
conductivity, and A := (k + l)/(2/c — 2) its contrast. In the circular setting, N coincident 
sources/receivers are evenly spaced on the circle of radius R and centered at the origin 
O between the angular range (0,7]. In the full- view case, 7 = 2tt while 7 < 2ir in the 
limited-view configuration. The position of s-th source (and r-th receiver) is denoted by x s 
(and x r , respectively) for s,r = 1 . . . N, with 6 S = js/N the angular position. We require 
that the circle is large enough to include the inclusion (R > 6). In the following, we set 
p := R/5 > 1. 

2.1 CGPTs and the linear system 

In the presence of D, the electrical potential u s resulting from a source at x s is given as the 
solution to the following conductivity problem [2]: 

V-((1 + (k-1) X d)Vu s )(x) = 0, xel 2 , 
u s (x) — T(x — x s ) = 0(|x| _1 ), \x\ — > +00, 

where T(x) = (l/27r) log|x| is the fundamental solution of the Laplacian in M 2 : AT(x) = 
So(x), with 5o being the Dirac mass at 0. 

Using asymptotic expansion of the fundamental solution, the MSR data V = (V sr ) Sjr 
being defined as V sr = u s (x r ) — F(x r — x s ), is linearly related to the GPTs of B as [4, 6]: 

* § \*\+\P\ 

Vsr= —^rd a T(z-x s )M a ^(X,B)d l3 T{z-Xr) + E sr + W sr , (2) 

where K denotes the highest order of GPTs in the expansion, E = {E sr ) S)T the truncation 
error (non-zero if K < 00), and W = {W S r)s,r the measurement noise following indepen- 
dently the same normal distribution: W sr ~ Af(0,a^ oisc ), of mean zero and variance o"n 0ise - 
The contracted GPTs, being defined as a harmonic combination of the GPTs [8], allow 
us to put (2) into an equivalent form [2]: 

Vsr = T (cosme s ,smmOs) (™T ( cosn9 Q r ) -J— +Esr + W sr , 

2irmp m v SJ \M% n M%J \smn9 r J 2^np n 

A sm M mn (A rn ) T 

(3) 
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where T denotes the transpose and the CGPT matrix M = (M- mn )m,n 1 has dimension 
2K x 2K. 

Recall that A = CD, with C being a N x 2K matrix constructed from the block 
C rm = (cos m9 r , sin m8 r ) and D a 2K x 2K diagonal matrix: 



/ Cn C12 
C21 C22 

\Cjvi Cn2 



CnkJ 



;D 



1 

2^ 



fh/p 



\ 



\ 



h/(2p 2 



l2/(Kp K )J 



(4) 



Here, I2 is the 2x2 identity matrix. With these notations in hand, we introduce the linear 
operator 



L(M) = CDMDC T , 

and rewrite (3) as: 

V = L(M) + E + W. 

In order to reconstruct M, we solve the least-squares problem: 

min ||L(M) - Vlli, 

M 



(5) 



(6) 



(7) 



where \\-\\f denotes the Probenius norm. It is well known that (7) admits a unique minimal 
norm solution M est = L^(V), with being the pseudo-inverse of L provided by the 
following lemma: 

Lemma 2.1. Let A,B be two real matrices of arbitrary dimension, and define the linear 
operator L(X) = AXB T . // A^,B^ are the pseudo-inverse of A,B respectively, then the 
pseudo-inverse of L is given by 



L\Y) = A t Y(B t ^ 



(8) 



Proof. This is a straightforward verification of the definition of pseudo-inverse, namely: 1) 
L^L and LL^ are self-adjoint; 2) LL^L = L and L^LL' = L'. For the first point: 

Lt(L(X)) = AtAX(BtB) T , 

which is self-adjoint since the matrices A^ A and B^B are symmetric by definition of pseudo- 
inverse; while for the second point, it follows from the definition again that 



L(L f (L(X))) = AA t AX(BB t B) T = AXB T = L(X). 
Similarly, one can verify the self-adjointness of LL^ and L^LL^ = iA 



□ 



1 Throughout the paper, we will write M mn for the m, n-th 2x2 building block, and (M) a j, for the a, 6-th 
entry in M. 
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2.2 Full-view setting 

In [2], we have investigated the resolving order of CGPT reconstruction in the full angle of 
view setting: 7 = 2tt. Given N > 2K, it has been shown that the matrix C is orthogonal 
(up to the factor N/2): 



2 

and the pseudo-inverse solution takes the form: 

Lt(V) = ^D^^VCD 1 . (9) 

Furthermore, the reconstruction problem is exponentially ill-posed. More precisely, the 
following result holds. 

Proposition 2.2. Let e a b be the 2Kx 2K matrix whose elements are all zero but the (a, b)th 
element is equal to 1. In the circular and full-view setting with N > 2K, the (a, b)-th singular 
value of the operator L, for a, b = 1, ... , 2K, is 

\ ab = N/(^ 2 \a/2-\ \b/2-\p^ 2 ^ b W ), (10) 

with the matrix e ab as the right singular vector, and i a f, = A~ fe 1 L(e a b) as the left singular 
vector. In particular, the condition number of the operator L is K 2 p 2 ^ K ~ 1 ^ . 

Proof. Using the fact that C T C = ^-1, we have, for any square matrices U and V, 

(L(U), L(V)> = ^ (DUD, DVD) , (11) 

where (•, •) is the termwise inner product. Since D is diagonal and invertible, we conclude 
that the canonical basis {e a b} a ^ is the singular vector of L, and the associated singular 
value is ||L(e a6 )|| F = ||De ab D|| F iV/2 = N/(8ir 2 \a/2] \b/2]p^ 2 ^ b ^). □ 

As a simple consequence, we have Lt(W) ab = A^W, / o6 ). When K is sufficiently 
large, the truncation error E is 0{p~ K ~ 2 ) and can be neglected if compared to W [2], and 
then by the property of white noise 



E(((M«% - (M) afe ) 2 ) < ^/E((Lt(W&)) = A-V noise) 

which is the result already established in [2]. Hence, it follows from (10) that the recon- 
struction of high order CGPTs is an ill-posed problem. Nonetheless the system has the 
remarkable property that low order CGPTs are not affected by the error caused by the 
instability of higher orders as the following proposition shows. 

Proposition 2.3. Let Mjf denote the CGPTs of order up to K , and let Jjk be the cor- 
responding linear operator in (3). Then, for any order K\ < K2 < N/2, the submatrix of 
L/j^ 2 (V) formed by the first 2K\ columns and rows is identical to the minimal norm solution 

Lk(V). 
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Proof. Let the N x 2K matrix Jk be the row concatenation of the 2K x 2K identity matrix 
I2K and a zero matrix. We have JJ^Jk = I2K and J^ i L^ 2 (V)Jx 1 is the submatrix of 
L^ 2 (V) formed by the first 2K\ columns and rows. Let and Ck be the matrices 
defined in (4). Because of (9), we have 

One can easily see that 
Thus, we have 

□ 

Numerically, can be implemented through either the formula (9) or the Conjugated 
Gradient (CG) method using (7). Simulations in [2] confirm that in typical situations, say, 
with K = 5 and 10% noise, the reconstructed CGPT is sufficiently accurate for the task 
such as the target identification in a dictionary. In the next section we present a location 
and orientation tracking algorithm for a mobile target based on the concept of CGPTs. 

3 Tracking of a mobile target 

At the instant t > 0, we denote by zt = [xt,yt] T £ the location and 6t £ [0, 27r) the 
orientation of a target Dt. 

D t = z t + R 6t D, (12) 

where R 9t is the rotation by t . Let M f be the CGPT of D t , and M D be the CGPT of D. 
Then the equation (6) becomes: 

V t = L(M 4 ) + Et + W t , (13) 

where Ef is the truncation error, and W; the measurement noise at time t. 

The objective of tracking is to estimate the target's location z% and orientation 0t from 
the MSR data stream Vf. We emphasize that these informations are contained in the first 
two orders CGPTs as shown in the previous paper [2]. Precisely, let Ax t = xt — Xt-i, 
Ayt = yt — Ut-x and AOt = 0t — Ot—ij then the following relations (when it is well defined) 
exist between the CGPT of D t and A-i [2]: 

NgCAVNftpt) = 2(Ax t + iAy t ) + e^N^A-ij/N^A-i), 

(14) 

Ng^MiCA) = 2(A^ + iAyt) + e^'N^CA-OMiCA-i). 

Hence when the linear system (14) is solvable, one can estimate zt,9t by solving and ac- 
cumulating Axt, Ayt and AOt- However, such an algorithm will propagate the error over 
time, since the noise presented in data is not properly taken into account here. 

In the following we develop a CGPT-based tracking algorithm using the Extended 
Kalman Filter, which handles correctly the noise. We recall first the definition of com- 
plex CGPT, with which a simple relation between Mj and M£> can be established. 
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3.1 Time relationship between CGPTs 

Let u = {l,i) T - The complex CGPTs NW,N( 2 ) are defined by 

= (AC - M™ n ) + i(M™ n + M% n ) = u T M mn u, 
= (M« n + ACJ + i(AC n - M-J = /M^u, 

where denotes the Hermitian transpose. Therefore, we have 

N (1) = U T MU and N (2) = U H MU, (15) 

where the matrix U of dimension 2K x K over the complex fields is defined by 



U 



(u ... 0\ 
u ... 

\0 ... uj 



(16) 



It is worth mentioning that and N^ 2 ) are complex matrices of dimension K x K. 

To recover the CGPT M mn from the complex CGPTs nW.N^, we simply use the 
relations 

1 1 

jurcc _S»nv(l) _l iv(2) \ yrcs _ _Ck(-]vr(l) , iy(2) \ 

? ? (17) 

M sc = -^CnW _ ]\K 2 ) ~) M ss = -WNW - ) 

AVJ -mn 2 ^ mn A, rrwi/> xyj -mn \ mn ^-^mnli 

where K, 9 are the real and imaginary part of a complex number, respectively. For two 
targets Dt,D satisfying (12), the following relationships between their complex CGPT hold 
[2]: 

N«(A) = ^ T N (1) p)F t , (18a) 
N^(D t ) = F t H N^(D)F t , (18b) 

where Ft is a upper triangle matrix with the (m, n)-th entry given by 

(F t ) mn =( n \xt + iyt) n - m e im9 K (19) 
\m J 

Linear operator T$: Now one can find explicitly a linear operator T< (the underlying 
scalar field is M) which depends only on zt,9t, such that Mt = T^Md), and the equation 
(13) becomes 

V t = L(T t (M D )) + E t + W t . (20) 

For doing so, we set J< := VFt, where U is given by (16). Then, a straightforward compu- 
tation using (15), (17), and (18) shows that 

M CC (A) = KJ f T M D KJ t , M cs (D t ) = &J?M D SJ U 
M sc (D t ) = %jjM D $lJ t , M ss {D t ) = 3?J t T M D 3?J t , 
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where M cc (D t ), M cs (D t ), M sc (D t ), M ss (D t ) are defined in (3). Therefore, we get the oper- 
ator T t : 

T t (M D ) = RU{RjjM D RJ t )R\J T + 3W($L/j T M D 3 J t )QU T + 

9U(SJ t T M D KJ t )KU T + 9U(9J t T M D 9J t )9U T = Mi . (22) 

3.2 Tracking by the Extended Kalman Filter 

The EKF is a generalization of the KF to nonlinear dynamical systems. Unlike KF which is 
an optimal estimator for linear systems with Gaussian noise, EKF is no longer optimal, but 
it remains robust with respect to noise and computationally inexpensive, therefore is well 
suited for real-time applications such as tracking. We establish here the system state and 
the observation equations which are fundamental to EKF, and refer readers to Appendix B 
for its algorithmic details. 



3.2.1 System state observation equations 

We assume that the position of the target is subjected to an external driving force that has 
the form of a white noise. In other words the velocity (V(r)) rg K+ of the target is given in 
terms of a two-dimensional Brownian motion (W / a (r)) rg j R + and its position (Z(r)) TgK + is 
given in terms of the integral of this Brownian motion: 

V(t) = V + a a W a (r), Z(r) = Z + f V{s)ds. 

Jo 

The orientation (©(T)) rgR + of the target is subjected to random fluctuations and its angular 
velocity is given in terms of an independent white noise, so that the orientation is given in 
terms of a one-dimensional Brownian motion (W / e(r)) rg R+: 

@(T) = e + a e We(r). 

We observe the target at discrete times tAr, t S N, with time step At. We denote zt = 
Z(tAr), vt = V(tAr), and 9t = ©(tAr). They obey the recursive relations 



a a (W a (tAr)-W a ((t-l)AT)), 

tAr 

W a {s) - W a ((t - l)Ar)ds, (23) 

(*-1)At 

c t = a e {W 9 (tAr) - W 9 {{t - l)Ar)) . 



at 
b t 



On 



v t -i + a t , 

Zt-i + vt-iAr + b t , 
t = 6 t -i + c t , 

Since the increments of the Brownian motions are independent from each other, the vectors 
(Ut)t>l given by 

U t 



are independent and identically distributed with the multivariate normal distribution with 
mean zero and covariance matrix 5] given by 



At 



/ a 2 a I 2 ^Arl 2 o\ 
4AtI 2 ^Ar 2 I 2 



V 



o 



(24) 
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The evolution of the state vector 




takes the form 



X t = FI W + U t , 



F = 



( 



I 2 o 0' 
Arl 2 h 
1 



(25) 



The observation made at time t is the MSR matrix given by (20), where the system 
state Xt is implicitly included in the operator T t . We suppose that the truncation error E t 
is small compared to the measurement noise so that it can be dropped in (20), and that 
the Gaussian white noise Wj of different time are mutually independent. We emphasize 
that the velocity vector v t of the target does not contribute to (20), which can be seen from 
(12). To highlight the dependence upon zt,9t, we introduce a function h which is nonlinear 
in zt, 0t, and takes as a parameter, such that 



Then together with (25) we get the following system state and observation equations: 



Note that (27a) is linear, so in order to apply EKF on (27), we only need to linearize (27b), 
or in other words, to calculate the partial derivatives of h with respect to Xt,yt,&t- 

3.2.2 Linearization of the observation equation 

Clearly, the operator L contains only the information concerning the acquisition system and 
does not depend on xt,yt,0 t . So by (26), we have 



h(X t ;M D ) = h(z t ,e t ;M D ) = L(T t (M D )). 



(26) 



X t = FX t _! + U u 

V t = h(X t ;M D )+W t . 



(27a) 
(27b) 




+ 



(28) 



(29) 



a a ,(KJ t T M D KJi) = ^{d Xt Jj)M D ^J t + mjM D m(d xt j t ) 
9 a ,(KJ f T M D 9Ji) = U(d xt jJ)M D %J t + RJ?M D Z(d Xt J t ) 

d xt {^jJw D m t ) = ^{d Xt jJ)m D m t + ^jJm D ^{d xt j t ) 
d Xt (Qj?M D Qj t ) = Q{d Xt jJ)M D Qj t + Qj^M D Q(d xt j t ) 



and d Xt Jt = Ud Xt Ft- The (m, n)-th entry of the matrix d Xt Ft is given by 




(30) 



The derivatives d yt T t (M.£)) and dg t T t (M.£)) are calculated in the same way. 
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3.3 Numerical experiments of tracking in the full-view setting 

Here we show the performance of EKF in a full angle of view setting with the shape 'A' as 
target D, which has diameter 10 and is centered at the origin. The path (zt,6t) is simulated 
according to the model (23) during a period of 10 seconds (At = 0.01), with parameters 
a a = 2, ao = 0.5, and the initial state Xq = (vq, zq, 6q) t = (— 1, 1, 5, — 5, 37r/2) T . We make 
sure that the target is always included inside the measurement circle on which N = 20 
sources/receivers are fixed, see Fig. 1. The data stream V t is generated by first calculating 
the MSR matrix corresponding to each Dt,t > then adding a white noise. 

Suppose that the CGPT of D is correctly determined (for instance, by identifying the 
target in a dictionary [2]). Then we use the first two orders CGPT M/j of D in (27b), and 
take (0, 0, 10, —0.5, 0) T as initial guess of Xq for EKF. 

We add 10% and 20% of noise to data, and show the results of tracking in Fig. 2 (a) (c) 
and (e). We see that EKF can find the true system state, despite of the poor initial guess, 
and the tracking precision decays as the measurement noise level gets higher. The same 
experiment with small target (of same shape) of diameter 1 is repeated in Fig. 2 (b) (d) and 
(f), where the tracking of position remains correct, on the contrary, that of orientation fails 
when the noise level is high. Such a result is in accordance with physical intuitions. In fact, 
the position of a small target can be easily localized in the far field, while its orientation 
can be correctly determined only in the near field. 



40 



30 
20 
10 


-10 



-20 
-30 
-40 

-40 -30 -20 -10 10 20 30 40 

Figure 1: Trajectory of the letter 'A' and the estimation by EKF. The initial position is 
(5, —5) while the initial guess given to EKF is (10, —0.5). The crosses indicate the position 
of sources/receivers, while the circle and the triangle indicate the starting and the final 
position of the target, respectively. In blue is the true trajectory and in red the estimated 
one. 
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Figure 2: Results of tracking using the configuration of Fig. 1 at different noise levels. First 
row: coordinate in x-axis. Second row: coordinate in y-axis. Last row: orientation. In the 
first column the target has size 10, while in the second column the target has size 1. The 
solid line always indicates the true system state. 
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4 CGPT reconstruction and tracking problem in the limited- 
view setting 

In this section we study the stability of CGPTs reconstruction and tracking problem in 
the case < 7 < 2ir, always under the condition that N > 2K, i.e., the number of 
sources/receivers is two times larger than the highest order of CGPTs to be reconstructed. 
Unlike in the full-view case, here C is no longer orthogonal in general, nonetheless one can 
still establish the SVD of L similarly as in Proposition 2.2. 

Proposition 4.1. Consider the concentric and limited-view setting with N > 2K , and 
suppose that C is of maximal rank. Let {fj. n } be the n-th largest eigenvalue of the matrix 
DC CD and let {v n } be the associated orthonormal eigenvector. Then the (a, b)-th singular 
value of the operator L is X a b = y/Jt^JIb, with the associated left singular vector the matrix 
Sab = v a v J ■ I n Particular, the condition number of the operator L is 

cond(L) = cond(DC T CD) < cond(C) 2 A- 2 p 2(i ^" 1) , (31) 
with cond(C) being the condition number of the matrix C. 
Proof. We first note that for any matrices U, V we have: 

(L(U),L(V)) = (U, (DC T CD)V(DC T CD)). 
Taking g afe = v a vj , and g a /y = v a >vj,, we get 

(L(g a& ),L(g a / 6 /)) = Ha'(v a vJ ,u a '4(DC T CD)) = fi a/ fi b >(v a vj ,v a/ vj,} 

= Saa'hb'/J-a/J-b, 

where 5 aa / is the Kronecker's symbol, which implies that ||L(g a {,)||^ = y^I^fe is the (a, b)-th 
singular value of L. We denote by /0 max (-), p m in( - ) the maximal and the minimal singular 
values of a matrix, then 

p max (DC T CD) = p max (CD) 2 < p m ax(C) 2 / o max (D) 2 , 
p min (DC T CD) = p min (CD) 2 > Anin (C) 2 p min (D) 2 , 

and the condition number of L is therefore bounded by cond(C) 2 ET 2 /? 2 ^ _1 \ □ 
4.1 Injectivity of C 

We denote by Vk the vector space of functions of the form 

K 

f(9) = ^e lk6 , (32) 

k=-K 

with Cfc £ C, and the subspace of Vk such that Co = 0. Functions of V K can be written 

as 

K 

f(e) = J2 a kcos(ke) + /3 k sin(ke), (33) 

k=l 
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with afc, (3k S C Observe that taking discrete samples of (33) at 9 S = js/N is nothing but 
applying the matrix C on a coefficient vector (a±, fi\ . . . ax, Pk)- We have the following 
result. 

Proposition 4.2. For any N > 2K , the matrix C is of maximal rank. 

Proof. Multiplying / G V K in (32) by e iKe ', and using the fact that cq = 0, we have 

K-\ 2K 
k=0 k=K+l 

K-l 2K-1 2K~l 

Ck-KC +2^ c k+i-i<e — °k e > («3 4 J 

fc=0 k=K k=0 

where c k = c^-k f° r k = 0, . . . , K — 1, and = e ie, Cfc + i_x for = K, . . . , 2K — 1. The 
iV vectors u s := (e )fe=o...2.ftr-i are linearly independent since they are the first 2K < N 
rows of a N x N Vandermonde matrix. Therefore, f{6 s ) = for s = 1 . . . N implies that 
dk = for all k = 0, . . . , — 1, which means that C is of maximal rank. □ 

Consequently, for arbitrary range < 7 < 2ir, a sufficient condition to uniquely deter- 
mine the CGPTs of order up to K is to have N > 2K sources/receivers. 

4.2 Explicit left inverse of C 

We denote by Dk(0) the Dirichlet kernel of order K: 

n rm - - ^({K + 1/2)9) 

D «^-X. K e - sin{e/2 ) ■ ( 35 ) 

We state without proof the following well known result about Vk- 

Lemma 4.3. The functions {Dk{0 — 2 2 r+i )}n=o,...,2K is an orthogonal basis ofVx- For 
any f,g€ Vk, the following identity holds: 



v 2K + 1 J y I 2K + 1 J ' 

n=l v 7 x ' 

where * denotes the complex conjugate. In particular, we have for n = 0, . . . , 2K 



(36) 



Lemma 4.4. Given a set of N > 2_fT different points < 6\ < . . . < 6^ < 2tt, there exist 
interpolation kernels h s 6 ^[jv/2J / or s = 1 ... TV, smc/i i/ia£: 

A 7 " 

f(P) = J2f(0s)h.(P) foranyfeV K . (38) 

s=l 
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Proof. When the number of points N is odd, it is well known [27] that h s takes the form 

A' 



t=l,t^s k 

When N is even, by a result established in [22] 



N 

sin 



sin 



(^) 



It is easy to see that in both cases /i s belongs to Vjjv/2j ■ d 

Now we can find explicitly a left inverse for C. 

Proposition 4.5. Under the same condition as in Lemma 4-4> we denote by h s the inter- 
polation kernel and define the matrix C = (C ks ) k , s as 

1 r2n 1 /"27T 

C 2 fe-i, s = -/ h s (9)cos(k9)d9, C 2k>s = - h s (9)sm(k9)d9. (41) 
TTien CC = I. In particular, if N is odd, the matrix C can 6e calculated as 

C»-i,. = ^ E fc. J cos ^_ j , C« t . = - g fc. (— J sm (-^ J • (42) 

Proof. Given u = (ai, /3i . . . olki^k) £ C 2 ^, and / the associated function defined by (33), 
we have (Cv) n = f{6 n ) for n = 1, . . . , N . Using (38) and (41), we find that 

»2tt 



(CC«) 2fc -i = - r f(9) cos(k9)d9 = a k , (43) 
(CC^) 2fc = - / W S m(k9)d9 = /3 k , (44) 

Jo 



and therefore, CCv = v. Observe that h s (9), cos(&0), and sin(k9) all belong to Vjjv/2J> so 
when iV is odd, we easily deduce (42) using (36). □ 

Remark 4.1. In general, the left inverse C in (41) is not the pseudo-inverse of C, and by 
definition, we have = C if CC is symmetric. If P v o(h s ) is the orthogonal projection of 
h n onto Vk, i-e., 

K 

C 2 fc-i, s cos(k9) + C 2fc;S sin(fc0), (45) 

fc=i 

then, P v o(h s )(9t) = (CC) S (. Therefore, C is the pseudo-inverse of C if and only if the 
interpolation kernel h s satisfies: 

P v o(h s )(9t) = P v o(h t )(e s ), for s,t = 1. . .N. (46) 
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Remark 4.2. Proposition 4.5 can be used in the noiseless limited- view case to reconstruct 
the CGPT matrix M from the MSR measurements V. In fact, from (5) it immediately 
follows that 

M = D 1 CVC T D _1 . 

This shows that in the noiseless case, the limited-view aspect has no effect on the recon- 
struction of the GPTs, and consequently on the location and orientation tracking. In the 
presence of noise, the effect, as will be shown in the next subsection, is dramatic. A small 
amount of measurement noise significantly changes the performance of our algorithm unless 
the arrays of receivers and transmitters offer a directional diversity, see Fig. 6. 

4.3 Ill-posedness in the limited-view setting 

We undertake a numerical study to illustrate the ill-posedness of the linear system (6) in the 
case of limited- view data. Fig. 3 shows the distribution of eigenvalues of the matrix C T C 
and DC T CD at different values of 7 with N = 101 and K = 50. In Fig. 4, we calculate the 
condition number of C T C and L (which is equal to that of DC T CD by (31)) for different 
orders K. From these results, we see clearly the effect of the limited-view aspect. First, the 
tail of tiny eigenvalues in Fig. 3. (a) suggests that the matrix C T C is numerically singular, 
despite the fact that C is of maximal rank. Secondly, both C T C and L rapidly become 
extremely ill-conditioned as K increases, so the maximum resolving order of CGPTs is very 
limited. Furthermore, this limit is intrinsic to the angle of view and cannot be improved by 
increasing the number of source/receivers, see Fig. 4 (c) and (d). 




10 20 30 40 50 60 10 20 30 40 50 60 

Order of CGPT Order of CGPT 



(a) Eigenvalues of C T C (b) Eigenvalues of DC T CD 

Figure 3: Distribution of eigenvalues (in log scale) of the matrix C T C (a) and DC CD 
(b). N = 101 sources are equally spaced between [0,7) on a circle of radius p = 1.2, and 
K = 50. Each curve corresponds to a different value of 7. The matrix C T C and DC T CD 
are calculated from these parameters and their eigenvalues are sorted in decreasing order. 
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Figure 4: Condition numbers (in log scale) of the matrix C T C (a) and the operator L (b) 
for different orders K between [1,50]. As in Fig. 3, N = 101 sources are equally spaced 
between [0,7) on a circle of radius p = 1.2. Fig.(c) and (d) are the same experiment as 
Fig. (a) and (b) but with N = 1001. 
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4.4 Reconstruction of CGPTs 



The analysis above suggests that the least-squares problem (7) is not adapted to the CGPT 
reconstruction in a limited-view setting. Actually, the truncation error or the noise of 
measurement will be amplified by the tiny singular values of L, and yields extremely instable 
reconstruction of high-order CGPTs, e.g., K > 2. Instead, we, in order to reconstruct 
CGPTs from the MSR data, use Thikhonov regularization and propose to solve 

min ]|L(M)-V]|i + u]|M]||, (47) 

M 

with /i>0a small regularization constant. It is well known that the effect of the regular- 
ization term is to truncate those singular values of L smaller than which consequently 
stabilizes the solution. The optimal choice of fx depends on the noise level, and here we 
determine it from the range [lO^KT 1 ] by comparing the solution of (47) with the true 
CGPTs. 

Here we reconstruct the CGPTs of an ellipse with the parameter N = 101, if = 50, 
and 7 varying between and 2tt. The major and minor axis of the ellipse are 1 and 0.5 
respectively. In Fig. 5 we show the error of the first 2 order CGPTs reconstructed through 
(47) and (7) at three different noise levels. It can be seen that, for small 7, the error 
obtained by (47) is substantially smaller. 




Figure 5: Error of reconstructed CGPT of an ellipse compared with true CGPT values at 
different noise levels. We solve (47) and (7) with N = 101, if = 50, and compare the first 
two orders with the true CGPT. The x-axis is the angle of view 7. Fig. (a): results of (47), 
Fig.(b): results of (7). 



4.5 Tracking in the limited-view setting 

The performance of the tracking algorithm can also be affected by the limited angle of 
view. We repeat the experiment of subsection 3.3 with 5 = 10, 7 = tt, and the same 
initial guess. In the first configuration, N = 21 sources/receivers are equally distributed 
between [0,7), see Fig. 6 (a). The results of tracking by EKF presented in Fig. 7 (a), (c) 
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and (e) show large deviations in the estimation of position, and a totally wrong estimation 
of orientation. In the second configuration, we divide the sources/receivers into 5 groups 
placed in a nonuniform way on [0, 2n), and each group covers only an angle range of 0.27T, 
see Fig. 6 (b). Although the total angular coverages are the same in both configurations, 
the second one gives much better tracking results, as shown in Fig. 7 (b), (d) and (f). 
These results clearly demonstrates the importance of a large angle of view (or a directional 
diversity) for the tracking problem. 



20- . 




- True trajectory _ 50 
Estimation 



-40 -30 -20 -10 



10 20 30 40 



- True trajectory 
Estimation 




-50 -40 -30 -20 -10 10 20 30 40 50 



Figure 6: Same experiment as in Fig. 1, with a limited angle of view 7 = tt. In Fig. (a) 
sources/receivers are equally distributed between [0,7), while in Fig.(b) they are divided 
into 5 groups. 



5 Conclusion 

In this paper we have provided a location and orientation tracking of a mobile target from 
MSR measurements in the full- and limited- view settings. Our algorithm is based on the 
concept of GPTs. In the limited-view case, the effect of noise is severe on the tracking. 
However, if the arrays of receivers and transmitters offer a good directional diversity, then 
satisfactory results can be obtained. It would be interesting to generalize our algorithms 
for tracking multiple targets. As a first step, a matching pursuit algorithm [21] would be 
appropriate for recognizing the targets. This will be the subject of a forthcoming work. 



A Kalman Filter 

The KF is a recursive method that uses a stream of noisy observations to produce an 
optimal estimator of the underlying system state [20] . Consider the following time-discrete 
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Figure 7: Results of tracking using the configuration of Fig. 6 at different noise levels. First 
row: coordinate in x-axis. Second row: coordinate in y-axis. Last row: orientation. First 
and second column correspond to the configuration in Fig. 6 (a) and (b), respectively. 
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dynamical system (t > 1): 

X t = F t X t -i + W t , (48) 
Y t = H t X t + V t . (49) 

where 

• Xt is the vector of system state; 

• Yj is the vector of observation; 

• Ft is the state transition matrix which is applied to the previous state Xt-i; 

• Ht is the observation matrix which yields the (noise free) observation from a system 
state Xt; 

• Wt ~ Af(0,Qt) is the process noise and Vt ~ Af(0,Rt) is the observation noise, with 
respectively Qt and Rt the covariance matrix. These two noises are independent 
between them, further, Wt of different time instant are also mutually independent 
(the same for Vt). 

Suppose that Xq is Gaussian. Then it follows that the process (Xt,Yt)t>o is Gaussian. 
The objective is to estimate the system state X t from the accumulated observations Y\-t := 
\Yi...Y t ]. 

The optimal estimator (in the least-squares sense) of the system state Xt given the 
observations Y\ : t is the conditional expectation 

x tl t = nX t \Y 1 ..t}. (50) 

Since the joint vector (X t ,Yi ]t ) is Gaussian, the conditional expectation IS Si linear 
combination of Y± : t, which can be written in terms of Xf_iu_i and Yt only. The purpose of 
the KF is to calculate x t u from x t -\\t-\ and Y t . 

We summarize the algorithm in the following. 

Initialization: 



Prediction: 



x |o=E[X ], P | = cov(Xq). (51) 

%-i = F t x t -i\ t -i, (52) 

Y t = Y t - #t%-i, (53) 

P t \t-i=FtP t -i\t-lFl + Qf (54) 



Update: 



S t = H t P t \ t _ 1 H? + R t , (55) 

K t = P t \t-iH?St\ (56) 

xt\t = xt\t-i + K t Y t , (57) 

Pt\ t = {I-K t Ht)Pt^\t-i- (58) 
To apply the KF algorithm the covariance matrices Qt,Rt must be known. 
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B Extended Kalman Filter 

Consider now a nonlinear dynamical system: 

X t = ft(X t - 1 ,Wt), (59) 
Y t = h t (X t ,V t ), (60) 

where Xt,Yt, Wt, Vt are the same as in the KF, while the functions ft, ht are nonlinear and 
differ entiable. Nothing can be said in general on the conditional distribution Xf|Yi : f due 
to the nonlinearity. The EKF calculates an approximation of the conditional expectation 
(50) by an appropriate linearization of the state transition and observation models, which 
makes the general scheme of KF still applicable [26]. However, the resulting algorithm is 
no more optimal in the least-squares sense due to the approximation. 

Let Fx = dx f (x t _iu i,Q), Fyy = dw f { x t-i\t-u 0)> the partial derivatives of / (with 
respect to the system state and the process noise) evaluated at (x t -i\t-i, 0), and let Hx = 
dxh(x t \ t _i, 0), Hy = dvh(x t \ t _i,0) be the partial derivatives of h (with respect to the 
system state and the observation noise) evaluated at (x t \t-i,0). The EKF algorithm is 
summarized below. 

Initialization: 

x |o=E[X ], P 0{0 = cov(X ). (61) 

Prediction: 

x t\t-l 
% 
Pt\t-i 

Update: 

St = HxPt\t-\H x + HvRtHy, 

% = %-i + K t Y t , 
P t \t = (I-K t H x )P^ Mt _ 1 . 
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